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m ; Abstract 

We present a simple analytical tool which gives an approximate insight into the 
stationary behavior of nonlinear systems undergoing the influence of a weak and 
\^ , rapid noise from one dominating source, e.g. the kinetic equations describing a ge- 

Q I netic switch with the concentration of one substrate fluctuating around a constant 

mean. The proposed method allows for predicting the asymmetric response of the 
genetic switch to noise, arising from the noise-induced shift of stationary states. The 
method has been tested on an example model of the lac operon regulatory network: a 
reduced Yildirim-Mackey model with fluctuating extracellular lactose concentration. 
^ . We calculate analytically the shift of the system's stationary states in the presence 

OO I of noise. The results of the analytical calculation are in excellent agreement with the 

I~^ ' results of numerical simulation of the noisy system. The simulation results suggest 

T^-j- ■ that the structure of the kinetics of the underlying biochemical reactions protects 

_il I the bistability of the lactose utilization mechanism from environmental fluctuations. 

O ' We show that, in the consequence of the noise-induced shift of stationary states, 

o 



cr: 



the presence of fluctuations stabilizes the behavior of the system in a selective way: 
Although the extrinsic noise facilitates, to some extent, switching off the lactose 
metabolism, the same noise prevents it from switching on. 
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1 Introduction 



Relatively simple biochemical systems regulated at the level of gene expression 
are capable of a complex dynamic behavior due to their intrinsic nonlinearity. 
The nonlinear kinetics of the biochemical regulation may result in various pat- 
terns of behavior, among which bistability is an extremely interesting one as a 
source of a switch-like behavior, a common strategy used by biochemical and 
cellular systems to turn a graded signal into an all-or-nothing response. An- 
other important feature associated with the bistability is hysteresis: in order 
to switch the system from one steady state to another, the input signal must 
surpass a given threshold. To switch back, the input signal must be decreased 
below another (smaller) threshold. This permits a discontinuous evolution 
of the system along different possible pathw ays, which may provide the sys- 
tem with an epigenetic (nongenetic) memory ( Laurent and Kellershohnl . Il999l ; 
Casadesus and D'Aril . 120021 : iFerrell . l2002h . 



Recently, growing attenti on has been focused on the study of stochastic as- 



pects of gene regulation (lAustin et al.l. l2006l: lElowitz et al.l. l2002l: iPaulsson 
200,4 iRosenfeld et~^ . l200,4 ISwain et al.l . I2OO2I : iTsimring et all . I2OO6I ). Fki^- 



tuations in a gene network are generally divided into 'intrinsic' and 'extrinsic'. 
This distinction depends on the point of view: we consider the low copy num- 
ber fluctuations in a single reaction (or set of reactions) under study as the 
'intrinsic noise', whereas the 'extrinsic noise' is connected with all the remain- 
ing, external processes which are not taken into consideration in detail. The 
extrinsic noise can originate from low copy number fluctuations in the reac- 
tions that are external with respect to the set of processes studied, as well as 
from other stochastically varying, unknown factors affecting our system. 



Reliable functioning of a cell may, on the one hand, require genetic networks 



to suppress or to be robust to fluctuations (ITabaka et al.l . l2008l : lElowitz et al. 



2002 



Becskei and Serranol . I2OOOI : lAlon et al.l . Il999l ). On the o ther hand, noise 



offers the opportunity to generate a switch-like behavior (IQzbudak et al 



2004 ) and a long-term heterogeneity in a clonal population (lElowitz et al 



2OO2I ). The presence of fluctuations in nonlinear systems such as genetic net- 
works may induce spontaneous switching between stationary st ates, emergence 

of ne w stationary states and disappearance of the existing ones ([Horsthemke and Lefever 
1984h . 



In this paper we present a simple analytical tool which gives an approximate 
insight into the stationary behavior of systems undergoing the influence of a 
weak and rapid noise from one nonlinear source, e.g. kinetic equations where 
the fluctuations of the concentration of one particular substrate dominate, 
and where that concentration enters into the equations in a nonlinear func- 
tion. The proposed method of mean noise expansion allows for predicting the 



noise-induced shift of stationary states which, in case of bistable systems, 
causes the asymmetric response to fluctuations. This asymmetry may, for ex- 
ample, facilitate switching off a genetic switch but prevent it from switching 
on. The occurrence of such an effect can give rise to the question of a possible 
optimization of the genetic switch for functioning in a noisy environment. 

We have tested our method on an example model of the lac operon. The lactose 
regulation system in the Escherichia coli bacteria is one of the most exten- 
sively studied examples of a biological switch: it allows for the maintenance of 
differences in the phenotype despite th e absence of genetic and environmental 
differences. Monod and Pappenheimer (JMonod et al.l . Il952l ) discovered the ef- 
fect of population heterogen eity on the level of an enti re bacterial population, 
whereas Novick and Weiner ( JNovick and Weineii . 119571 ) identifled the same ef- 
fect on the level of individual E. coli cells: in the same external conditions, the 
bacteria were either able or unable to metabolize lactose. The studies of the 
expression of /3— galactosidase (the enzyme which breaks down the lact ose into 
a simpler sugar) (iNovick and Weineii . 119571 : ICohn and Horibatal.ll959l) and of 



the d irect lac gene transcription activity at the cellular level (jOzbudak et al. 



20041 ) show that the cells can be in one of two discrete states: either fully 
induced, with the transcription (and, consequently, enzyme) levels reaching 
a maximum, or uninduced, with negligible transcription and enzyme levels. 
The induction may be triggered by applying even a quite short stimulus: a 
temporary increase in the extracellular lactose level. 



Different nonlinear dynamical models of the underlying chemical kinetics were 
proposed to expl ain the origins of the switch-like behavior of the lacto se uti- 
lization network (iLaurent and Kellershohnl . Il999l : lOzbudak et al.l . 120041 ) . This 
direction of r esearch has led to a rnore de tailed model proposed by Yildirim 
and Mackey (lYildirim and Mackeyl . 120031 ). It explicitly incorporates all the 
relevant biochemical processes along with experimentally motivated kinetic 
constants, and, tested on empirical data, displays a good agreement with ex- 
periments. According to this model, the switch-like behavior of the lactose 
operon results from the bistability of the kinetic equations. 



Since the changes of the extracellular lactose concentration are the primary 
factor which controls the induction and uninduction of the lactose metabolism 
in E. coli, we focus our attention on this, completely external, process influenc- 
ing the lac operon system. Within the example model based on the Yildirim- 
Mackey framework, we analyze how the weak and rapid Gaussian fluctuations 
in the extracellular lactose concentration (and their different intensity) affect 
the lac gene expression. We do not take into account the intrinsic fluctua- 
tions (modeling their effects deserves a separate study) but our analysis may 
be the flrst step to the interpretation of the experimental measurements of 
stochasticity in lac operon expression, in terms of the discrimination between 
the effects of the intrinsic and extrinsic noises, which is itself a challenging 



task. It is worth noting that lElowitz et alJ (120021 ) have shown that in systems 



consisting of several reactions (in particular, also in the lac operon system in 
E. coli) the extrinsic noise often gives a much stronger contribution to the 
gene expression than the intrinsic fluctuations. 

Using the proposed method of mean noise expansion, we analytically calculate 
the noise-induced shift of the stationary states of the model, which gives rise to 
the asymmetric response of the system to fluctuations: the effective stabiliza- 
tion of the uninduced state and the destabilization of the induced state. We 
show that the results of the analytical calculation are in excellent agreement 
with the mean stationary states obtained from the numerical simulation of the 
noisy system. We also show the consequences of that shift: Varying the noise 
intensity, we measure mean times of the transition between the uninduced 
and induced states. In this way we check when the system becomes resistant 
to the fluctuations and when, on the contrary, the fluctuations facilitate the 
switching between those states. 

The paper is organized as follows: in Section 2 we present the analytical 
method of calculation of the noise-induced shift in stationary states of a sys- 
tem. In Section 3 the example model of the lac operon is described. Section 4 
presents the results: the application of the analytical method compared to 
the numerical results (Subsection 4.1), and changes in the mean times of the 
transition between the uninduced and induced states being the consequence 
of the noise-induced shift (Subsections 4.2 and 4.3). 



2 Theory 



The general method of treatment of dynamical systems undergoing the influ- 
ence of a weak and rapid noise from one nonlinear source, which we present be- 
low, can be applied, for example, to the models of genetic regulatory networks. 
A genetic switch should be described by the equations of chemical kinetics, 
i.e. neglecting all sources of noise, except one: a fluctuating concentration of a 
substrate which enters into the equations as a parameter. The concentration 
of that substrate should fluctuate weakly but rapidly around a constant mean. 



Assume that: 
a) The system is described by a set of stochastic differential equations: 

^ = F(X,MP,)). (1) 

X can denote here the vector of concentrations of reactants. Pt is a param- 
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eter, for example, a concentration of a substrate, which does not depend on 
X. 

b) Pt is a stochastic process, fluctuating in time t around the mean P and 
having a constant variance cr^ -C 1 (weak fluctuations). 

c) The fluctuations of Pt are rapid enough not to correlate with the time scales 
of the processes described by the Sys. (1). The characteristic time Tgys of the 
system is given by -m^, where A is the greatest eigenvalue of the Jacobian 
of (1) (a standard linearization procedure). The characteristic time scale of 
the process Pt is determined by its correlation time r. Therefore, r <^ r^ys- 

d) Pt enters into the system (1) only in the function h{Pt) only. 

e) The deterministic system 

^ = F(X,MP)) (2) 

with a constant parameter P equal to the mean of Pt has steady states 

X*(P). 

If b) and c) are fulfilled, we can assume that the trajectories of the stochastic 
system (1) fluctuate weakly and rapidly around a certain constant average 
(X(Pt)). This means that the behavior of the system is quasi-stationary, i.e. 
even if the probability density of X(Pj) has more than one maxima, the tran- 
sitions between them are very unlikely. Therefore we will consider only the 
trajectories which fluctuate around one of the maxima: (X(PJ) will be then 
the position of that maximum. Since the fluctuations are weak, the maxima 
of the probability density of X(P() are close to the steady states of the deter- 
ministic system (2). 

The response of the system to noise in the parameter P will be a shift of 
the mean, around which the trajectories of the stochastic system fluctuate, 
by a small value of A with respect to the corresponding steady states of the 
deterministic system: 

(X(P,)) = X*(P + A). (3) 



In order to flnd A, we take the noisy trajectories of (1) at a certain time t and 
average them over the ensamble of realizations: 

= (F(X(P,),MP,))) (4) 



At weak and rapid fluctuations, we can assume that X (indirectly depending 
on Pt) experiences only an averaged contribution from the noise, so it can be 
replaced by the constant (X(Pt)), which depends on the mean P, but also on 
other parameters of the process Pt. For example, if Pj is Gaussian, i.e. fully 
characterized by its mean and variance, (X(Pt)) = const{P,a'^) (it depends 



on the mean value of the parameter and on the strength of its fluctuation). 
Then, the only term in (4) which directly depends on the fluctuating values 
of Pt is h{Pt). The averaging in (4) can be thus separated: 

= F((X(P,)), (MP,))). (5) 



At a given time t, Pt = P + 6P for each trajectory. Assuming that the devi- 
ations from the mean are small, we make a Taylor expansion of h{Pt) around 
the mean P up to the second order: 

{h{Pt)) = {h{P + 6P)) = h{P) + h'{P){6P) + h"{P){6P') + ... (6) 



The average deviation from the mean (SP) = 0. The mean square deviation is 
equal to the variance: {SP"^) = o^ . (If Pt is a Gaussian process, then also the 
third-order term {6P^) = 0). The Eq. (6) is now approximated by: 

(MPt)) = HP) + lh"{P)a\ (7) 



The deterministic steady states in (3) are given by the equation: 

= F(X*(P + A), h{P + A)), (8) 

whereas (5) is now 

= F((X(P,)), hiP) + h"iP)a'). (9) 

Using (3), we can compare (8) and (9). The approximate value of the shift A 
can be therefore calculated from 

h(P + A) = h{P) + -h"{P)a\ (10) 

For the weak noise, the shift should be small, so it can be further approximated 
by the Taylor expansion of h{P + A) up to the flrst order in A. The Eq. (3) 
then takes the form: 

(X(P,)) = X*(P+^^a2), (11) 

iih'{P) 7^0. 
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Fig. 1. Schematic representation of the lactose operon regulatory system. 

The approximation lies in the fact that we attribute the shift of the stationary 
states to the functional dependence of h{Pt) only. In some cases, the multi- 
plicative dependence on the system's variables, h{Pt) f (X.{Pt)) , may be also 
important. Here we neglect it assuming that the noise is weak and rapid, so 
that the system variables can be treated as constant. In case our method does 
not show the noise-induced shift of stationary states (for example when h" = 0, 
i.e. the noisy parameter enters linearly into F), it may still turn out that a 
noise-induced shift is present, n amely because the multip l icativ e dependence 
of the noise becomes important ( JHorsthemke and Lefeven . Il984l ). 



3 Example: lactose svi^itch model 



Below we present an example model of lac operon regulatory network to which 
our method will be applied. 



3. 1 The regulation mechanism of the lac operon 



The scheme of the feedback regulation rn e chani sm of the induction in the lac 
operon in E. coli (JYildirim and Mackeyl . l2003l ) is presented in Fig. 1: The 
extracellular lactose (Lg) is transported through the cell membrane by the 
enzyme permease (P). The intracellular lactose (L) is then broken down into 
glucose, galactose, and allolactose {A) by the enzyme /3-galactosidase {B). By 



a positive feedback loop, the presence of allolactose enables the production of 
permease and /3-galactosidase enzymes. 

The lac operon consists of three genes, lacZ, lacY, and lacA preceded by 
a promoter/operator region. The lacZ gene encodes for the mRNA respon- 
sible for the production of /3-galactosidase, whereas the lacY gene produces 
the mRNA for permease. (The third gene, lacA, producing thiogalactoside 



transacetylase, does not play a role in the lac operon regulation (iBeckwith 



19871 ).) The transcription of this part of DNA is controlled by a repressor 
protein (R). If the inducer (allolactose) is absent, the repressor binds to the 
operator DNA sequence and makes the transcription of genes by the RNA 
polymerase impossible. In the presence of inducer, a complex is formed be- 
tween allolactose and the repressor that prevents the latter from binding to the 
operator region. The RNA polymerase is then able to initiate the transcription 
of the lacZ, lacY, and lacA genes to produce mRNA (M). Subsequently, the 
mRNA is translated into the apropriate enzymes (P and B). 

3.2 Yildirim-Mackey model 

Our study is based on the model of lac op eron regulation proposed by Yildirim 



and Mackey (j Yildirim and Mackeyl . 120031 ). The model consists of five equations 



of chemical kinetics for the reactions described in Subsection 3.1, involving 
mRNA (M), allolactose (A), lactose (L), /3-galactosidase (B) and permease 

(P): 



aM fi{e-^^^' Ar,,) + To - 1m M 



dM 
dt 

^ =aB e-^^^ M,^ - 7B P 

f = «A P g,{L) -PaB MA) -^aA (12) 

f = ax P h{Le) -PlP 92{L) -aAB g,{L) ~%L 

^ =ap e-'^(-^+-^) M.^+.^ - 7p P 



fi{A) describes the effector-repressor dynamics of the transcription enhanced 
by allolactose: 



The square term derives from the average number of allolactose molecules 



required to inactivate the repressor, which (according to lYildirim and Mackey 



(120031 ): lYagil and Yagill (Il97ll )) is approximately equal to 2. /2, gi, g2 and h 



are rates of the irreversible Michaelis-Menten reactions: 

/2(^) = K^^ 9i{L) = j^, 

a and [3 denote the gain and loss rates for the reactions. Ki is the equilibrium 
constant for the repressor-allolactose reaction. K2 is the equilibrium constant 
for the operator-repressor reaction, K = 1 + K2Rtot, and Rtot is the total 
amount of the repressor, r represents the delays associated with the finite 
time required to complete the transcription (tm) and the translation (rp and 
tb). For example, A^^ = A{t — tm)- The 7 = 7 + /! are the coefficients for 
the terms representing decay of species due to chemical degradation (7) and 
dilution (/x). The exponential factors take into account the dilution of mRNA 
due to cell growth. Even if allolactose is totally absent, on occasion repressor 
will transiently not be bound to the operator and RNA polymerase will initi- 
ate transcription. Although the mRNA production rate dM/dt would be then 
nonzero (a leakage transcription), it is necessary to add an empirical constant 
Fn to the model to obtain a l eakage rate that agrees with experimental values 



( JYildirim and Mackeyl . 120031 ). A more detailed derivation of the above equa- 



tions, an elaborate estimation of the parameter values, as well as the results of 



testing the model on experimental data can be found in (JYildirim and Mackeyl . 



20031 ). The model reproduces the bistable behavior of the lactose operon for 



realistic extracellular lactose concentrations. 



3.3 Reduced model 



In this paper, we study a reduced model derived from the full Yildirim-Mackey 
equations. The simplified model is less time-consuming for numerical simula- 
tions (Subs. 4.3 and 4.2) and easier for analytical treatment (Subs. 4.1) but, 
at the same time, it captures all the important features of the full model. 

Our model consists of three equations of kinetics for mRNA (M), allolactose 
(A) and lactose (L) concentrations in the bacterial cell: 

dM ^ l+Ki A^ I -p r, /\# 

-It -^M i+K^Rtot+K, A^ + -L - 7M iW 

^ =^BM(a^^-/3^^)-7^A (15) 

f =fcpM(ai^^-/?^^^)-a^A:5M^-7LL 
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Fig. 2. Steady states of the system (for mRNA concentration). Marked are the initial 
and final points of the simulations described in Section 4. 

The model has been obtained from Eqs. (12) by the assumption that: i) The 
enzyme concentration is proportional to the mRNA concentration. The as- 
sumption is based on the fact that the amount of proteins observed in procar- 
iotic cells is in general proportional to the amount of t heir transcripts (whic h 
is, in particular, the case for the lac operon in E. coli) ([Mathews et all 119991 ) . 
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(16) 



ii) The time delays are small enough to fullfi ll: e"^'^"' ^ 1, and iii) X(t + T ^) ~ 
X{t). The values of parameters estimated by (lYildirim and Mackeyl . l2003l ) (see 
Table A.l) allow for such assumptions. Assumption iii) is very well fullfilled in 
the neighborhood of stationary points, but also the differences between X{t) 
and X(t + Tx) in transient regions are not large compared to the distance 
between stationary points (see the trajectories in Figs. 4, 5). 



The stability properties of the model reproduce very well those of the ful l 
Yildirim-Mackey model (see Fig. 2 and compare with (lYildirim and Mackeyl . 
20031 )): for the concentrations of the extracellular lactose 



0.0277 mM < L, < 0.0618 mM, 



(17) 



the system has two steady states (induced and uninduced). Assuming that 
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the extracellular lactose concentration is constant and does not fluctuate in 
time, the model shown above is deterministic. If L^ lies in the range of bista- 
bility, then, after a sufficiently long time, the concentrations of allolactose, 
lactose and mRNA reach one of the steady states. Switching between the two 
steady states is then impossible unless the extracellular lactose concentration 
changes (increasing or decreasing Lg beyond the bistable region would enable 
a hysteretic response of the system) . 



3.J^ Fluctuations in the extracellular lactose concentration 



We investigate how the fluctuations in the extracellular lactose concentration 
induce the shift of stationary states of the system and spontaneous transi- 
tions between those states. To model the fluctuations in Le, we have chosen 
Gaussian noise, assuming that the natural environment, where wild-type bac- 
teria remain, undergoes the influence of very many random factors. If their 
contributions are independent, they add together into the noise of a Gaussian 
distribution around a certain mean. We therefore assume that the noise in L^ 
is Gaussian and fast with respect to othe r processe s in th e system. We model 



it using the Ornstein-Uhlenbeck process (JGardinen . l2004l ) 
dL 



dt 



-eiL,-L,) + ^m (18) 



which fluctuates in time around the mean value Lg. ^(t) is a Gaussian white 
noise of intensity 7 and autocorrelation (C(t)^(s)) = 6(t — s). The correlation 
time of the fluctuations in Lg, tqu = V^ = 2 ■ 10^^ min = 1.2 s, has been 
chosen significantly larger than the time step in the numerical calculations 
{6t = 2 ■ 10"'^ min) but smaller than the fastest time scale of the system 
Tsys ~ 10^^ min (see Appendix A). In this way we obtain a rapidly fiuctu- 
ating stochastic process which does not correlate with the time scales of the 
biochemical processes described by the model. Otherwise, if the noise were as 
slow as the system's reaction to it, the A, L, and M concentrations would have 
enough time to adapt to the instantaneous values of Lg. In that case, the fur- 
ther analytical treatment described in this paper would be impossible. Taking 
into account the niobility of the E. coli bacteria (mean velocity ~ 30/im/s, see 



e.g. iDiLuzio et al.l ( 120051 )). granularity of the intestinal content and motions 
of intestinal villi, we may suppose that the fiuctuation rapidity assumed here 
i s realistic . The variance of the extracellular lactose fiuctuations is 7^/(2^) 



(lGardineii . l2004l ). The value of 7 controlling the noise intensity will be varied 



in our simulations. 
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4 Results 



Our goal was to determine how the model of the lac gene expression is sensitive 
to fluctuations in the extracellular lactose concentration. Using the proposed 
method of noise expansion, we calculated analytically the noise-induced shift of 
the stationary states of the system (15) with the external lactose concentration 
Le perturbed by the Gaussian noise given by the Eq. (18). 

In order to show the effects of the asymmetric response of the system to 
noise, we performed numerical simulations of the trajectories of the system 
(see Appendix B for detailed technical information). We measured the mean 
time of noise-driven transitions between the induced and uninduced states for 
different values of the mean extracellular lactose concentration Le and with 
different noise intensities 7 (see Appendix B for the simulation details). To 
show the fluctuation strength on the graphs in a more intuitive way, the noise 
intensity unit used in the below-presented flgures is the standard deviation of 



Le, equal to ^77(2^). 

4.1 Analysis of noise-induced stabilization and destabilization 



To analytically predict the behavior of the steady states of the system under 
the influence of noise, we use the approximate method of mean noise expan- 
sion presented in Sec. 2. Treating the instantaneous value of Le(t) as a small 
perturbation from the mean value, we expand the noise in L^ around its mean 
and average the equations of kinetics over an ensemble of possible trajecto- 
ries, to get the mean stationary concentrations of A, L, M in the noisy system. 
The fluctuating parameter Le enters into the equations (15) in the Michaelis- 
Menten form h{Le) = ^ ^^^ . In the Fig. 3 we compare the noise-induced 
shifts of stationary states of the mRNA concentration: (i) obtained from the 
numerical simulations, (ii) calculated by numerical solution of the Eq. (10), 
(iii) and by the further approximation (11) which here takes the form of 

2/i'(Le) 9iKL^ + L,y ^■' 



All the results are in excellent agreement: positions of the stationary states 
in the bistable regime shift right with increasing noise intensity, which results 
in the destabilization of the induced steady state for small Le and the sta- 
bilization of the uninduced steady state for large L^. The steady-state mean 
mRNA concentrations for a given value of Le decrease due to noise (compare 
with Fig. 5). 



12 



o 

%— ' 

CD 

c 

0) 

o 

c 
o 
o 



E 

CO 

c 
g 

^— ' 
en 



10" 



10" 



1 ■ ■ ■ ■ ^^-^s^l" 


. . ■ ■ 1 


■ 


■ 


1 


■ 


^..-^'""''^ 


/>' 


3.5 


- 


+ 




^,-f^ 1^180 


// ■ 






1 




/J; CD 

f '^130 

1 ^ 




3.2 
2.9 


' y 




- 


■\ 80 


\V 1 1 


2.6 




1 












■ X^ 0.028 0.03 0.032 0.055 


0.06 


^\^^ Le[mM] Le 


[mlVl] 


- K/\*i\ ^ 


^--^....^^^^ 


- 


■ Ivl (Lg) 


: mean M, stoch. simulation + ^"""~""'*---.^ 




■ r/i /I , A\ 




---^, 


IVI ^Lg+Zi^ 


^y| /I Mn2-h"/lOh'\\ 




^> 


IVI ^Lg+CJ 11 /^^ll )) 




.,-H— t 


-— -r*^ 





0.03 



0.04 0.05 

Lg [mM] 



0.06 



Fig. 3. Stationary states {M{Le)) for mRNA concentration at different mean con- 
centrations of extracellular lactose, found using the approximate method of noise 
expansion. Under the influence of noise (o"^ = 0.01 mM), the stationary states shift 

right by A with respect to the deterministic stationary states M*{Le). The dashed 

2 
lines for A calculated from the Eq. (10) and for A = , ^ ^ ^ (Eq. (11)) overlap. 

Comparing with the results of the numerical simulations (+), one can note that the 
approximation is excellent. 

4.-2 Transition from the induced to uninduced state 



In order to show the consequences of the noise-induced shift of the stationary 
states, we performed a series of simulations in which the system started in 
the induced state. The initial state of the system was given by the determin- 
istic (i.e. calculated in the absence of noise) stationary concentrations of A, L 
and M at the induced state for a given mean extracellular lactose level L^ (see 
Fig. 2 and Table B.l). The results (Fig. 4) show that the system is generally re- 
sistant to the fluctuations in the extracellular lactose concentration. However, 
the noise slightly destabilizes the system, i.e. it enables switching from the in- 
duced to uninduced state. Such transitions are only possible for Le very close 
to the bifurcation point Lg = 0.0277 (the beginning of the bistable region). A 
small increase in the mean extracellular lactose concentration causes a drastic 
increase of the mean transition time. For example, after changing the mean 
extracellular lactose concentration from Lei = 0.0277 mM to Le2 = 0.0279 
mM, the noise intensity (measured by variance 7^/(26*)) must be almost one 
order of magnitude stronger to force the transition to occur in a similar time. 
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Within the simulation time (1000 min) no backward transitions have been 
recorded. The graph of mean transition times vs. noise intensity shows that 
the stronger the fluctuations are, the shorter is the mean transition time i.e. 
the easier is the uninduction. As we will see in the next subsection, the noise- 
driven induction depends on the noise strength in quite a difl^erent manner. 
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Fig. 4. Transition from the induced to uninduced state. Simulation of changes in 
mRNA concentration, depending on the strength of fluctuations in extracellular 
lactose concentration (here, for clarity, measured by the standard deviation of L^). 
Top panel: Example realizations (single trajectories) for mean extracellular lactose 
concentrations Lgi = 0.0277 mM (see Fig. 2) and different noise intensities. Bottom 
panel: Mean times of the transition for Lei = 0.0277, Ze2 = 0.0278 and Zgs = 0.0279, 
depending on the noise intensity. The initial and final points for the trajectories are 
same as in Fig. 2). A slight increase in the mean extracellular lactose concentration 
from -Lei to Le2 causes a drastic increase of the mean transition time. 
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Fig. 5. Transition from uninduced to induced state. Simulation of changes in mRNA 
concentration, depending on the strength of fluctuations in extracellular lactose 
concentration. Top panel: Example realizations (single trajectories) for mean ex- 
tracellular lactose concentration Le5 = 0.0620 (beyond the bistable region) and 
different noise intensities. We have chosen the starting points same as the coordi- 
nates of the steady state at the bifurcation point Le4 = 0.0618 (see Fig. 2). For 
stddev(Le) = 0.010, the fluctuations make the induction impossible (also, the mean 
mRNA concentration shifts down). In other cases, the noise causes a delay of the 
induction. Bottom panel: Mean times of the transition beyond the bistable region 
for Les = 0.0620, L^q = 0.0630 and Lg? = 0.0650, depending on the noise intensity. 
We have chosen the starting points same as the coordinates of the steady state at 
the bifurcation point, and the final points as the deterministic steady states (see 
Fig. 2). Here, the stronger the fluctuations are, the longer is the mean transition 
time, i.e. the more difficult is the induction. 
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4-3 Transition from the uninduced to induced state 



Next, we performed simulations which started from the uninduced state (See 
Fig. 2 and Table B.l). Here we observe that, on the contrary to the previous 
simulations, the fluctuations of the extracellular lactose concentration stabi- 
lize the uninduced state (see Fig. 5). In the bistable region (increasing Le up 
to the bifurcation point Le4 = 0.0618, i.e. to the end of the bistable region), 
even very strong fluctuations cannot force the induction. Beyond the bistable 
region (Les = 0.0620, Lgg = 0.0630, Lg? = 0.0650), we have chosen the start- 
ing points same as the coordinates of the steady state at Le4, and flnal points 
as the deterministic induced steady states for Les, Le6 and Ley. In this case, 
the noise causes a delay of induction, and sufficiently strong fluctuations make 
the induction impossible at all, whereas the mean mRNA concentration even 
decreases. The dependence of the mean transition time on the noise intensity 
is opposite than in the case of uninduction: here, the stronger the fluctua- 
tions are, the longer is the mean transition time, i.e. the more difficult is the 
induction. 



5 Summary 



We present a simple analytical tool which gives an approximate insight into 
the stationary behavior of nonlinear systems undergoing the influence of a 
weak and rapid noise from one dominating source. These can be, for example, 
the kinetic equations describing a genetic switch with the concentration of 
one substrate fluctuating around a constant mean. The method is applicable 
in those cases when the noise enters into the equations in a form of a nonlin- 
ear function of a parameter which fluctuates weakly and rapidly, so that the 
fluctuations can be assumed to contribute to the state of the system as an 
average (i.e. when the noise is rapid enough not to directly couple with the 
system's dynamics). The proposed method of mean noise expansion allows for 
predicting the effect of the asymmetric response to noise, which, for example, 
may facilitate switching off a genetic switch but prevent it from switching on. 

The method has been tested on the example mo del of the lac operon regula - 



tory network: a reduced Yildirim-Mackey model (JYildirim and Mackeyl . 120031 ) 
subject to fluctuations in extracellular lactose concentration Lg. We analyzed 
how the fluctuations, modelled by the Gaussian noise, and their different inten- 
sities affect the lac gene expression predicted by the model. Using the method 
of mean noise expansion, we calculated analytically the shift of the system's 
stationary states in the presence of the fluctuations. We have shown that the 
results of the analytical calculation are in excellent agreement with the mean 
stationary states obtained from the numerical simulation of the noisy system. 
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The shift of stationary states gives rise to the asymmetric response of the 
system to the noise: the effective stabilization of the uninduced state and the 
destabilization of the induced state. We show the consequences of that shift by 
analyzing the mean times of the numerically simulated noise-driven transitions 
between induced and uninduced states in the bistable region (i.e. in the range 
of extracellular lactose concentration Lg where both induced and uninduced 
states are possible). The simulation results show that the system as a whole is 
highly resistant to fluctuations. It can be stochastically induced only for the 
lowest possible Lg, i.e. very close to the left boundary of the bistable region, 
and only at a high noise intensity. On the other hand, stochastic induction is 
impossible in the bistable region even for the highest possible Lg, i.e. in the 
closest neighborhood of the right boundary of the bistable region. Moreover, 
noise inhibits induction even beyond the bistable region: if fluctuations in Lg 
are strong enough, the lactose metabolism remains uninduced, although at 
weaker fluctuations or in their absence it would switch to the induced state. 
Thus, in the result of the steady-state shift, the influence of noise on the 
induction and uninduction of lactose metabolism is quite opposite: whereas 
the spontaneous uninduction (although difficult) is facilitated by noise, the 
same noise causes the suppression of the induction. 

The above results suggest that the bistability of the lactose utilization mecha- 
nism is "protected" by the structure of the kinetics of the underlying biochemi- 
cal reactions: switching the lactose switch due to an extracellular fluctuation is 
very difficult. Moreover, the "protection level" is different for the induction and 
uninduction: the possibility of random switching on the lactose metabolism 
is much more strongly protected than the possibility of random switching it 
off. Although, on the mathematical level, such a behavior results from the 
structure of the kinetic equations, a question may be posed if this particu- 
lar shape of kinetics can have a deeper explanation, for example whether the 
"protection" of the switch against the external fluctuations is connected with 
preventing an unnecessary energetic effort? This suggests a new direction of 
study: an analysis of the noisy system of lactose metabolism in E. coli from 
the energetic point of view. 

It is worth noting that the noise-induced shift of stationary states and the 
emergence of new stationary states in systems with multiplicative noise (such 
as the system under study) are the effects which may be described in terms 
of the stationary probability distributions. These can be obtained by solv- 
ing the Fokker-Planck equation associated to the n oisv equations of kinetics , 



treat e d as a multi-dim ensional Langevin equation ( IHorsthemke and Lefeven . 



1984 : iGardineii . l2004l ) . Further study in this direction seems to be an inter- 
esting perspective. The effect of intrinsic fluctuations as well as the response 
of the Yildirim-Mackey model to slow extrinsic noise also deserves a separate 
study. 
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A Details of the reduced Yildirim-Mackey model 



Table A.l presents the parameters of the model (same as used in (JYildirim and Mackey 
20031 )) and the data in the Table A. 2 shows the stationary states of the system 
for two example values of external lactose concentration, close to bifurcation 
points. Table A. 2 



To 


7.25 X 10"^ mM/min 


/i 


0.0226 min"^ 




dA 


1.76 X 10^ min-^ 


tb 


2.0 min 




as 


1.66 X 10^2 min"^ 


tm 


0.1 min 




ai 


2.88 X 10^ min"^ 


rp 


0.83 min 




au 


9.97 X 10"^ mM/min 


K 


7.2 X 10^ 




ap 


10.0 min"^ 


Ki 


2.52 X 10^ mM" 


-2 


Pa 


2.15 X 10^ min-1 


Ka 


1.95 mM 




/3l 


2.65 X 10^ min~^ 


Kl 


0.97 mM 




lA 


0.52 min"^ 


Kl. 


0.26 mM 




IB 


8.33 X 10"^ min-i 


Kl, 


1.81 mM 




IL 


0.0 min-^ 


kB 


0.677 




iM 


0.411 min"^ 


kp 


13.94 




IP 


0.65 min"^ 









Table A.l 

Parameters of the model (same as used in |Yildirim and Mackeyl . l2003l )) 







Stationary states [mM]: 






Uninduced 


Induced 


Le[mM] 


A* 


L* M* 


A* 


L* M* 


0.0278 


4.00 • 10-3 


9.40 • 10-2 2.12 • 10-6 


1.04-10-2 


0.129 7.52 - 10-5 


0.0610 


1.22 • 10-2 


0.216 3.19-10-6 


0.386 


0.280 7.90 - 10-^ 



Table A. 2 

Stationary states of the system for two example values of external lactose concen- 
tration, close to the bifurcation points. Table A. 3 presents the characteristic times 
of the system calculated for the stationary points given in Table A. 2. 



20 



Characteristic times [mini 



Lp = 0.0278mM 



Lp = O.OeiOmM 



Uninduced 



Induced 



Uninduced 



Induced 



n T2 



•^3 



n T2 



T3 



n T2 



T3 



Tl 



T2 T3 



1.4 3.1 13 



43 0.64 0.42 



1.1 6.4 41 



0.21 2.7 0.055 



Table A. 3 

Characteristic times of the system calculated for the stationary points given in Ta- 
ble A.2 



B Simulation details 



We performed the numerical simulations of the trajectories of the system de- 
scribed by the Eqs. (15) with the lactose concentration Lg perturbed by the 
Ornstein-Uhlenbeck noise given by t he Eg. (18) . The equations were s olved 
numerically using the Euler scheme (jPress et al.l . Il993l : iMannellal . 120021 ) with 
the timestep 6t = 2 ■ 10~^. The simulations were run starting from one of 
the deterministic steady states for a given Lg. The time was measured until 
the stochastic system's trajectory gets into a close neighborhood (of a radius 



D = \/bA? + bl? + (5M2 = 0.005 mM) of the other deterministic stationary 
state. The initial and final points of the simulations are presented in Table B.l. 
The simulation was monitored to prevent fluctuating concentrations become 
negative but no such event occured during the simulations. The number of 
simulation runs was A^ = 1000. 



Le[mM] 


Ai\mM\ 


Li[mM] 


Mi[10-5mM] 


Af{mM\ 


L/[mM] 


M/[10-5mM] 


Lei = 0.0277 


0.0969 


0.128 


7.521 


0.00398 


0.0937 


0.212 


Ze2 = 0.0278 


0.104 


0.129 


8.600 


0.00400 


0.0940 


0.212 


Ze3 = 0.0279 


0.109 


0.129 


9.406 


0.00400 


0.0946 


0.212 


Ze4 = 0.0618 


0.0141 


0.224 


0.360 


- 


- 


- 


Ze5 = 0.0620 


0.0141 


0.224 


0.360 


0.393 


0.285 


80.8 


Ze6 = 0.0630 


0.0141 


0.224 


0.360 


0.399 


0.289 


82.5 


Ze7 = 0.0650 


0.0141 


0.224 


0.360 


0.412 


0.298 


85.9 



Table B.l 

Initial (i subscript) and final (/ subscript) concentrations of allolactose (^), lactose 
(L), and mRNA (M) for simulations described in Section 4. 
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